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The method of averaging is used to investigate the phenomenon of capture into resonance for 
a model that describes a Keplerian binary system influenced by radiation damping and external 
normally incident periodic gravitational radiation. The dynamical evolution of the binary orbit while 
trapped in resonance is elucidated using the second order partially averaged system. This method 
provides a theoretical framework that can be used to explain the main evolutionary dynamics of a 
physical system that has been trapped in resonance. 



1. INTRODUCTION 



In two previous papers |]l|j2| , we considered the long term nonlinear evolution of a Keplerian binary system that is 
perturbed by a normally incident periodic gravitational wave, and in a recent work Q we considered the additional 
effect of radiation damping, which is of interest in connection with the observed behavior of the Hulse- Taylor binary 
pulsar PSR 1913+16 These studies have been concerned with the issue of gravitational ionization, i.e. the 

possibility that an external periodic gravitational wave could ionize a Keplerian binary system over a long period 
of time. The impetus for this subject has come from the close conceptual analogy between gravitational ionization 
and a fundamental physical process, namely, the electromagnetic ionization of a hydrogen atom. That is, in these 
studies one hopes to learn about the disposition of gravitational radiation to transfer energy and angular momentum 
in its interaction with matter. In our recent investigation we encountered an interesting dynamical phenomenon 
connected with the passage of the binary orbit through resonance. As the binary system loses energy and angular 
momentum by emitting gravitational waves, its semimajor axis and eccentricity decrease monotonically on the average; 
however, this process of inward spiraling could stop if the system is captured into a resonance. The resonance 
condition fixes the semimajor axis; therefore, if the semimajor axis decreases to the resonant value and the orbit is 
trapped, it then maintains this value on the average while the external perturbation deposits energy into the system 
to compensate for the radiation damping on the average. It turns out that along with this energy deposition, the 
external tidal perturbation can also deposit angular momentum into the binary orbit so that its eccentricity decreases 
considerably during the passage through resonance. This was the situation for the particular instance of resonance 
trapping reported in [pi . In general, the orbital angular momentum can increase or decrease while the orbit is trapped 
in resonance. Figure y] depicts passage through a (1 : 1) resonance in which the eccentricity decreases. A similar 
phenomenon but with an increase in eccentricity has been reported in the recent numerical study of the three-body 
problem by Melita and Woolfson ||. The same situation can occur over a long time-scale in a (4 : 1) resonance in 
our model as depicted in Figure |[ The dynamical phenomena associated with an orbit trapped in a resonance occur 
over many Keplerian periods of the osculating ellipse; therefore, it is natural to average the dynamical equations over 
the orbital period at resonance. This partial averaging removes the "fast" motion and allows us to see more clearly 
the "slow" motion during trapping. It is possible to provide a description of the slow motion as well as a theoretical 
justification for the transfer of angular momentum by means of the method of second order partial averaging near a 
resonance. The purpose of the present paper is to study this phenomenon theoretically using the method of averaging 
for the resonances that occur in a Keplerian binary that is perturbed by the emission and absorption of gravitational 
radiation. 

Let us consider the simplest model involving a perturbed Keplerian (i.e. nonrelativistic) binary that contains the 
effects of radiation reaction damping and external tidal perturbation in the lowest (i.e. quadrupole) order. The tidal 
perturbation could be due to external masses and gravitational radiation; in fact, we choose in what follows a normally 
incident periodic gravitational wave for the sake of simplicity The Keplerian binary orbit under the combined 
perturbations due to the emission and absorption of gravitational radiation in the quadrupole approximation is given 
in our model by the equation of relative motion 
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where x(t) is the relative two-body orbit, r = |x|, and k — Go(mi + mj). The binary consists of two point masses 
77ii and 77i2, n%i + 777,2 = M, that move according to X\(t) = (777,2 /M) x(i) and Xa(i) = — (mi/M) x(i) , so that the 
center of mass of the system is at the origin of the coordinates. In fact, the center of mass of the binary can be taken 
to be at rest in the approximation under consideration here. The explicit expressions for TZ and K, are 
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where v = x(£), an overdot denotes differentiation with respect to time, a and (3 are of the order of unity and are the 
amplitudes of the two independent states of linear polarization of the normally incident wave, p is a constant phase, 
and fl is the frequency of the external wave. 

It is interesting to transform the dynamical system (|l])-([|) to dimensionless form. To accomplish this, let x — » i?ox 
and t — > Tot where Rq and To are scale parameters. Under this transformation, k — > kTfi / Rq and K. remains unchanged 
if we let fl — > CI/Tq. Let us further restrict i?o an d T by the relation fcT 2 = i?Q, so that we can set k = 1 in the 
dynamical equations; for instance, this condition is satisfied if Rq is the initial semimajor axis of the unperturbed 
Keplerian orbit and 2ttTq is its period. Furthermore, the dynamical system (Q)-(|3|) is planar; therefore, it is convenient 
to express these dimensionless equations in polar coordinates (r, 9) in the orbital plane. The result is 
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where 5, < 5 « 1, is the dimensionless strength of radiation reaction and is given by 
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while e, < e << 1, is the dimensionless strength of the external periodic perturbation. In this paper, we let S — eA, 
where A,0<A<oo, is a parameter that is fixed in the system; in this way, we avoid dealing with a two parameter 
(e, S) perturbation problem. In particular, we consider only perturbations that correspond to fixed directions from 
the origin of this parameter space. The full two parameter problem would require the consideration of perturbations 
corresponding to all curves in the parameter space. 

In the absence of radiative perturbations (e = and 5 = 01, the dynamical system (|]) describes a Keplerian ellipse. 
It is therefore useful to express the dynamical equations (Q) in terms of Delaunay variables that are closely related 
to the elements of the osculating ellipse. This is the ellipse that the relative orbit would describe at time t, if the 
perturbations were turned off at t. The osculating ellipse always has the same focus, which is taken to be the origin of 
the (polar) coordinates in the space of relative coordinates. Let the state of relative motion be described by (x, v), or 
equivalently (r, 0,p r ,pg), at time t; then, the energy of the motion fixes the semimajor axis a of the osculating ellipse 
and its eccentricity is subsequently fixed by the orbital angular momentum pg . Only two angles are left to determine 
the osculating ellipse completely: the orientation of the ellipse in the orbital plane given by g and the position on the 
ellipse measured from the periastron given by the true anomaly v. The latter is obtained from p r pe — esint) and the 
equation of the ellipse 
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— = 1 + ecosu, 



and the former from 9 — v. The relevant Delaunay "action-angle" variables (L,G,£,g) are thus defined by 



L:=a 1/2 , G:=pe = L{\ 



e 
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£ := u — esinw, g:=0 — v, (6) 

where u is the eccentric anomaly of the osculating ellipse, r — a(l — e cos it), and £ is the mean anomaly. The dynamical 
system (|j) in terms of Delaunay variables is given briefly in Appendix [X] and used in the following section. 

The Delaunay equations of motion are useful for the investigation of periodic orbits using the Poincare surface of 
section technique It has been shown in that nearly resonant periodic orbits exist in system (0) for sufficiently 
small e and A. These correspond to (to : 1) resonances, where (to : n) refers to the resonance condition muj = nQ. 
Here to and n are two relatively prime integers and lo = 1/L 3 is the Keplerian frequency of the orbit. A linear 
perturbation treatment fliO[ first revealed resonant absorption at (to : 1) resonances. There could, in principle, be 



other periodic orbits whose existence is not revealed by our method 1 1 1 

In our numerical investigation of the simple nonlinear model described above, we found |3j instances of resonance 
trapping during which the behavior of the osculating orbit could not be inferred in a simple manner on the basis of 
equation (Q). However, the dynamics of the averaged equations in resonance is simpler to analyze and it turns out 
that our numerical results |3| can be adequately explained using the second order partially averaged dynamics. The 
phenomenon of resonance trapping appears to be of basic significance for the origin of the solar system; therefore, it 
is worthwhile to develop a general theoretical framework for the study of the evolutionary dynamics while trapped in 
resonance. 

The dynamics of a system when it is locked in resonance is interesting in any circumstance involving more than one 
degree of freedom; for instance, let us suppose that the resonance condition fixes an action variable — say, energy — and 
for a one dimensional motion this would then imply that the state of the system at resonance is definite. However, 
if other action variables are present, they will not necessarily remain fixed while the system is trapped in resonance. 
Instead, the state of the system will in general vary and its dynamics at resonance is best investigated using the 
method of averaging. This is a generalization of the simple procedure that is commonly employed in Hamiltonian 
dynamics: The Hamiltonian is averaged over certain "fast" variables and the resulting averaged Hamiltonian is used 
to derive new dynamical equations that presumably describe the "slow" motion in a certain averaged sense. The 
general method is described in Appendix [B|, and it is applied to the dissipative dynamical system under consideration 
here in the rest of this paper. 

Resonance is a general and significant physical phenomenon and the description of the state of a physical system 
while trapped in resonance is of intrinsic importance. The inherent dynamics at resonance is trivial for a one di- 
mensional oscillator, but is rich in physical consequences for higher dimensional systems. While the general methods 
described here could be applied to a wide variety of physical problems, we confine our attention to a single model. Our 
results may, however, be of qualitative interest in dealing with the three-body problems that arise in the discussions 
of the origin of the structure in the solar system. 

The present paper relies on the results of our recent work ||. We have repeated here only what is needed for the 
discussion of the dynamics at resonance; for a more precise and complete presentation of the background material, 
our papers should be consulted. 

Finally, a basic limitation of our model should be noted. The only damping mechanism that we take into account 
is gravitational radiation reaction; this is consistent with our assumption that the binary consists of Newtonian point 
masses moving in vacuum except for the presence of background gravitational radiation. In this model, a theorem of 
Neishtadt can be used to show that resonance trapping is a rather rare phenomenon However, taking due account 
of the finite size and structure of astronomical bodies and the existence of an ambient medium, we would have to 
include in our model — among other things — the additional damping effects of tidal friction as well as the various drag 
effects of the ambient medium and electromagnetic radiation (cf., for instance, fl|l2|-^5|). These additional frictional 
effects could well combine to violate the condition N of Neishtadt's Theorem P,M. Thus, resonance trapping may 
not be so rare in astrophysics after all [O] . Inclusion of these additional effects is beyond the scope of our work. 



2. PARTIAL AVERAGING NEAR A RESONANCE 

We will consider the dynamics of the model system (Q) that is derived in ||. It has the following form when 
expressed in the Delaunay elements for the Kepler problem under consideration (cf. Appendix |A|): 
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L = -e—+eAf L , 

G=-e^+eAf G 
dg 

t- .^ + .A/„ (7) 
where eH ey ± is the Hamiltonian corresponding to the external perturbation and 

n cxt = ^ 2 [aC{L, G, £, g) cos fit + /3S{L, G, £, g) cos{flt + p)] . (8) 

Here 

5 °° 
C(L, G,l,g) = 2° 2e2 cos2 3 + a 2 ^(A, cos 2<? cos - 5„ sin 2g sin 2/£), 

5 °° 
<S(L, G,£,g) = -jO 2 ^ 2 sin2g + a 2 ^{A v sin2#coS2/£ + 5^ cos 2g sin i/f), 

Av = ~J^2 [ 2v < 1 ~ e ') J > e ) " ( 2 ~ e ') J -M] > 

B„ = -^(1 - e 2 ) 1 / 2 [eJ>e) - - e 2 )J v {ve)] , (9) 

J„ is the Bessel function of the first kind of order v, and e = (L 2 — G 2 ) 1 / 2 /L. The radiation reaction "forces" /d, 
D G {L,G, £,<?}, are certain complicated functions of the Delaunay variables given in Appendix [a|. In fact, we will 
only require the averages of the fr> given by 

Id := ^ f D (L,G,£,g)d£. 
27r Jo 

These have been computed in [||, and are given by 

fG = -j^(& + 7e 2 ), 

ft = 0, fg = 0. (10) 

In order to study the dynamics of the system (Q) at resonance, we will apply the method of averaging. We note here 
that averaging over the fast angle £ and the time t gives the correct approximate dynamics for most initial conditions 
via Neishtadt's Theorem as explained in our previous paper ||. Here we are interested in the orbits that are captured 
into resonance. To study them, we consider partial averaging at each resonance. 

Let us fix the value of L at the (m : n) resonance, say L — with mjL\ — nfl, and consider the deviation of L 
from the resonance manifold. To measure this deviation, we introduce a new variable T> given by 

L = L^+e 1/2 V (11) 

and a new angular variable if by 

1 nfl 

£ = T t + if = — t + If. 

L* m 

The scale factor e 1 / 2 ensures that, after changing to the new variables in (Q), the resulting equations for V and ip have 
the same order in the new small parameter e 1 / 2 and therefore the system is in the correct form for averaging. These 
new coordinates are standard choices in the mathematical literature (for more details see, for example, |l6| ] or p8|). 
It is important to emphasize here that the small parameter in the actual dynamics is e; however, the small parameter 
turns out to be e 1 / 2 in this case for the averaged dynamics. 
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To effect the coordinate transformation, we use the expansion 



1 1 

and find 
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1 - e 1/2 — + e—j- + 0(e 3/2 ) 



(12) 



V = -e 1/2 Fu - eVF 12 + 0(e 3 / 2 ), 
G = -eF 22 +0(e 3 / 2 ), 



g = eF 42 + 0(e 3 / 2 ), (13) 

where the F^ (G, nVlt/m + tp, g, t) are defined such that the first index refers to the equation in which it appears and 
the second index refers to the perturbation order in powers of e 1 / 2 that is employed. These quantities are given by 

F U := ^(L*, G, —t + tp, g, t) Af L (L*,G, —t + tp, g), 
at m m 

Fa := -j- j£ (L., G, —t + tp, g, t) - A^ L*, G, —t + tp, g), 

F 22 := — 2^(i„ G, —t + tp, g, t) - Af G (L*,G, —t + <p,g), 
og m m 

F32 := £*, G, — t + 5 , i + A/ £ L„ G, — t + tp, «? , 

^42 := ^(L.,G, —t + tp, g, t) + A/ a (L*,G, — * + v , 5 ). (14) 

The system dl3| ) is 27rm/f2 periodic in the temporal variable — since Ti. cxt is 27r/f2 periodic in time — and is in 
time-periodic standard form. Anticipating our intention to average to second order, we will apply an averaging 
transformation (for a detailed exposition see |18|). It is the characteristic property of this transformation that it 
automatically renders system (O) in a form such that to lowest order the new system is exactly the first order 
averaged system and the second order averaged system can be simply obtained by averaging the new system (cf. 
Appendix [b]). To obtain the desired transformation, we define 

Fij := / (G, — s + tp, g, s) ds, 

Z7TTO Jo m 

and the deviation from the mean for F±i 

n n 

\(G,<p,g,t) :=F 11 (G,—t + p,g,t)~F 11 . (15) 
m 

Furthermore, we define A(G, tp, g, t) to be the antiderivative of A(G, tp, g, t) with respect to t with the additional 
property that the average of A should vanish, i.e. 

2-xvn/Q, 

A(G, tp, g, s) ds = 0. 

Moreover, we note that both A and OA/dtp have zero averages. Our averaging transformation is given by 

V = V-e 1 / 2 A(G,$,g,t), G = G, tp = tp, g = g. 

The averaging transformation is constructed so that its average becomes the identity transformation. 
Let us observe that if G, tp, and g depend on t as solutions of the system (|l3]), then 

A x i/2f 3V \ dA s 
A = A - £/ (lt)^ + 0(£) - 



After applying the averaging transformation, we find that the system (113) takes the form 
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G = -eF 22 + 0(e 3 / 2 ), 
?= eF 42 +0(e 3 / 2 ). 
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0(e 3 / 2 ), 



(16) 



Finally, we drop the 0(e 3 / 2 ) terms in d||) and average the remaining truncated system to obtain the second order 
partially averaged system 

T> = -e 1/2 F n -eVF 12 , 
G = —eF22, 



F 
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(17) 



This system is the averaged form of system ( p"3| ) after dropping its 0(e 3 / 2 ) terms; however, this coincidence is fortuitous 
in this case. In general, one has to employ an averaging transformation in order to obtain the second order averaged 
system. 

To explain the evolutionary dynamics at resonance, we will replace the actual dynamical equations by the second 
order partially averaged system. As explained in Appendix this is a reasonable approximation over a limited 
time-scale. We remark that although the second order partially averaged system will be used to explain some of the 
features of our model system near its resonances, the actual dynamics predicted by our model is certainly much more 
complex than the averaged equations reveal. In particular, we expect that near the resonances — perhaps in other 
regions of the phase space as well — there are chaotic invariant sets and, therefore, transient chaotic motions || . The 
averaged system is nonlinear and could in general exhibit chaos; however, we have not encountered chaos in the second 
order averaged equations obtained from the model under consideration here. 

It remains to compute pjy, where Fij are defined in (]l~4|). For this, we recall equation (||) and set £ — nVLt/m + (p, 
expand the trigonometric terms in the Fourier series using trigonometric identities, and then average over the variable 
t. The required averages involving the external perturbation can be computed from the average of 7^ C xt. In fact, after 
some computation, we find that 
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H oxt {L*,G, — t + tp,g,t)dt 



T C (L*,G, g) cos rrup + T S (L*, G, g) sin mip, 



where, for n = 1, 
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[ — aB m (e) sin 2g + (3A m (e) sin 2g sin p + (3B m (e) cos 2g cos p] , 



(18) 



and for 1, 7i C xt = so that in this case we can define T c = T s = 0. 

The averages of Jd, D G {L, G,£,g}, are given in (|Io|). The terms involving radiation damping that appear in the 
partially averaged system are obtained from these expressions by Taylor expansion about the resonant orbit using 
equation (til). For example, we will use 
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The second order partially averaged system (|17|) is thus given explicitly by 
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(19) 



where we have dropped the tildes. It is clear that L in the expressions involving T c , T s , and e must be replaced by 
L*, the value of L at resonance. 

Having derived the equations for the averaged dynamics (|l9|), we now turn our attention to the consequences of 
these equations and the comparison of predictions based on them with the actual dynamics given by (Q) . 



3. FIRST ORDER AVERAGED DYNAMICS 



The first order partially averaged system, obtained from (n9|) by dropping the O(e) terms, is given by 



V 



,1/2 



mT c sin rmp + mT s cos mip — AT(G) 



G = 0, 



9 = 



(20) 



In this approximation, the variables G and g are constants fixed by the initial conditions while the remaining system 
in V and ip is equivalent to a pendulum-type equation with torque; namely, 



(p H j(mT c sin my — mT s cos my) — j A T(G). 



(21) 



We also have a second order differential equation — a harmonic oscillator with slowly varying frequency — for the 
deviation T> given by 



where 



V + eufV = 0, 



r3m 2 m . . ] 1 / 2 

. (T c cos my + T s sm my) 
T* 



(22) 



(23) 



These results can be formally justified if they are recast in terms of a new temporal variable given be e 1 / 2 ^ however, 
we use t in order to facilitate comparison with the actual dynamics. 

To show that (^||) is an oscillator, we must show that to is a real number. To this end, we suppose that during 
capture into the resonance the orbit is near an elliptic region of the first order partially averaged system ( p0| ) . This 
system is Hamiltonian with an effective energy of the form 



,1/2 



where U represents the effective potential energy given by 

U := —(T c cos my + T s sin my) + (Ar)y. 
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If (£>, ip) — (0, (fo) is an elliptic rest point of the first order partially averaged system; then, U'(<po) = and U"(<po) > 0, 
where a prime denotes differentiation with respect to ip. It follows that 

U"((po) = m 2 (T c cos mipo + T s sinnupo), 

so that wo, w := [3U" '(^o)] 1 ^ 2 / '1%, must be real. 

To show that the frequency £ := e^-^w is slowly varying, we differentiate this frequency with respect to time to 
obtain 

9m 3 

£ = -e — — (T s cosrmp - T c sinrmp)V + 0(e 3/2 ). 
2L°w 

It follows from the first order averaged equations that tp — po is expected to be of order unity; thus, £ is nearly constant 
over a time-scale of order e" 1 / 2 since £ is proportional to e. In particular, X> varies on the time-scale e 1 ' 2 *. Inspection of 
equations and ( p2| ) reveals that while T> is predominantly a simple harmonic oscillator with frequency £o = e 1 ^ 2 ^, 
the motion of tp in time could be rather complicated involving essentially all harmonics of £o- Therefore, ip cannot in 
general be considered a slow variable in time. 

Librational motions (periodic motions in the phase plane) of the pendulum-type averaged equation correspond 
to orbits of the original system that are captured into the resonance. On the other hand, if there are no librational 
motions, then all orbits pass through the resonance. Thus, we observe a necessary condition for capture: the pendulum 
system must have rest points in its phase plane. The rest points of the first order averaged system are given by V = 
and if = ipo, with 

Rsm(nup + n) = -AY, (24) 

where mT c = Rcosr] and mT s = — Rsinn. As an immediate consequence, we have the following proposition: For 
the (m : n) resonance, if n ^ I, there is no capture into resonance. On the other hand, for n = 1, if there are 
librational motions, then we must have A|T/i? < 1. These observations suggest that in the presence of sufficiently 
strong radiation damping, i.e. for A sufficiently large, there are no librational motions and, as a result, each orbit will 
pass through all the resonant manifolds that it encounters. In particular, since L = a 1 / 2 , the semimajor axis of the 
corresponding osculating ellipse will collapse to zero. Moreover, capture is only possible for (m : n) resonances with 
n = 1 when A\T/R\ < 1. In this case, if A|r/i?| < 1, then the rest points in the phase plane of the pendulum-type 
equation (^f[) at the resonance are all nondegenerate. This is precisely Neishtadt's condition B. 

The quantities Y and R depend (nonlinearly) on the variables (L, G, g) as well as the parameters of the system; 
therefore, the precise range of their ratio Y/R is difficult to specify in general. However, the value of this ratio can 
be determined numerically. In fact, to find orbits that are captured into resonance as displayed in Figures 1-3, we 
use the first order averaged system. A rest point of the first order system corresponds to a "fixed" resonant orbit 
with T> = and constant G, g, and ip = ipo. The main characteristics of these orbits do not change with respect to 
the slow variable e 1 / 2 ^; that is, the resonant orbits are essentially fixed only over a time-scale of order e~ x / 2 . After 
choosing G and g, we solve for ipo at a rest point, and then test to see if the resulting resonant orbit corresponds to 
a libration point of the first order averaged system (e.g. Wq must be positive). If it does, we convert the point with 
coordinates (L», G, (fo, g) to polar coordinates and then numerically integrate our model (|]) backward in time from 
this starting value. When the backward integration in time is carried out over a sufficiently long time interval, the 
orbit is expected to leave the vicinity of the resonance. After this occurs, we integrate forward in time to obtain an 
orbit of (^) that is captured into the resonance. 

The first order partially averaged system (^0|) is Hamiltonian; therefore, we do not expect that its dynamics would 
give a complete picture of the average dynamics near a resonance for our dissipative system. In particular, the 
librational motions of the pendulum-like system are not structurally stable. Indeed, the phase space for the full four 
dimensional system (|2p| ) is foliated by invariant two dimensional subspaces parametrized by the variables G and g. If 
a librational motion exists in a leaf of the foliation, then the corresponding rest points, the elliptic rest point and the 
hyperbolic rest point (or points) at the boundary of the librational region for the associated pendulum-like system, 
are degenerate in the four dimensional phase space of the first order averaged system; their linearizations have two 
zero eigenvalues corresponding to the directions normal to the leaf. In the second order partially averaged system, 
viewed as a small perturbation of (p0|), we expect that only a finite number of the degenerate rest points survive. 
These correspond to the continuable periodic solutions described in |3j. We expect the corresponding perturbed rest 
points to be hyperbolic. As a result, they persist even when higher order effects are added. The dimensions and 
positions of the stable and unstable manifolds of these rest points determine the average motion near the resonance. 
In particular, if one of these rest points is a hyperbolic sink, then there is an open set of initial conditions that 
correspond to orbits permanently captured into resonance. Our numerical experiments have provided no evidence for 
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this behavior. However, another possibility — that is consistent with our numerical experiments — is that a rest point of 
the perturbed system in phase space is stable along two directions but unstable in the remaining two directions. Thus 
a trajectory with initial point near the corresponding stable manifold will be captured into resonance and undergo 
librational motions near the resonant manifold until it spirals outward along the unstable manifold. 



4. DYNAMICAL EVOLUTION DURING RESONANCE 



The dynamics of system ( |19[ ) is expected to give a close approximation of the near resonance behavior of our 
original model (Q) over a time-scale of order e -1 / 2 . In particular, a basic open problem is the following: Determine 
the positions of the rest points and the corresponding stable and unstable manifolds of ([[9]) . A solution of this problem 
would provide the information needed to analyze the most important features of the behavior of the orbits that pass 
through the resonance as well as those that are captured by the resonance. Unfortunately, rigorous analysis of the 
dynamics in the four dimensional phase space of system ( |l9"| ) seems at present to be very difficult. Thus, in lieu of a 
complete and rigorous analysis of the dynamics, we will show how to obtain useful information from approximations 
of the second order partially averaged system. 

A typical example of the behavior of the orbit as it passes through a (1 : 1) resonance is depicted in Fig. [l[ The 
semimajor axis undergoes librations of increasing amplitude around the resonant value. Furthermore, the eccentricity 
of the orbit generally decreases while the orbit is trapped. 

We propose to analyze the oscillations of V by using the 0(e) approximation of the second order differential 
equation for T> derived from the second order partially averaged equations (|T9|). In fact, a simple computation yields 
the following differential equation 

V + e 7 X> + ew 2 V = 0, (25) 

where 

7 := m {~^ C0Smi P ~ sinm^ + ^p^g (146 + 37e 2 ) (26) 

is evaluated at L = L r in equation (|25| ) and w is given by equation (p3|). Similarly, the second order partially 
averaged system ( |l9| ) can be used to derive an equation for the temporal evolution of ip to order e; however, the 
resulting equation turns out to be identical with equatio n (2l| ). 

It is clear from inspection of the differential equation (J25 ) that L = L* + t x l 2 T> oscillates about its resonant value 
L* with a libration frequency given approximately by £o = e 1 ' 2 Wo. The magnitude of this frequency is in agreement 
with our numerical results to within a few percent for the case of resonance trapping depicted in Fig. [TJ. Moreover, 
the amplitude of the oscillations will increase (decrease) if 7 < (7 > 0). 

The sign of 7 varies with the choice of parameters, variables, and the order of the resonance. Our numerical 
experiments — which have by no means been exhaustive — indicate that the (1:1) resonance for the linearly polarized 
incident wave with a = 1 and (3 = is special in that 7 is negative at the elliptic rest points of the first order averaged 
system. However, for a resonance with order m > 1, 7 is not of fixed sign. Nevertheless, our numerical experiments 
indicate that the amplitude of librations generally increases over a long time-scale (cf. Fig. ||). Figures || and || 
illustrate the apparent richness of the evolutionary dynamics near resonances with m > 1 . It should be emphasized 
that the damped (or anti-damped) oscillator ( p5| ) approximates the actual dynamics only over a time-scale of order 
e" 1 / 2 (cf. Appendix^). 

To obtain an approximate equation for the envelope of the oscillations, we consider a time-dependent change of 
variables for the oscillator (|25| ) defined by the relation 

V = V e - (e/2) J"o T(s)ds . (27) 
After a routine calculation, we obtain the differential equation 

V + ew 2 V = 0(e 3/2 ). (28) 



Thus, to order e, V(t) in equation (27) is the solution of a harmonic oscillator equation with a slowly varying frequency 
£ = e 1 / 2 ?!!. In fact, equation (E8f) is identical to the oscillator (22) 
averaged system. The solution V[t) of the second order equation (25) 
the amplitude of V. 



with frequency £ obtained from the first order 
in this approximation is obtained by modulating 



We note that if 7 and w are constants, then formula (E7J) reduces to 
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V{t) = Ae'^ t/2 cos(e 1/2 wt + r), 



where A and r are constants depending on the initial conditions; this is the standard result for the damped (or 
anti-damped)harmonic oscillator with constant coefficients. 

Equation fl27|) certainly agrees qualitatively with the numerical experiments reported in Figures To check the 
accuracy quantitatively, we numerically integrated system (^) and simultaneously evaluated the integral of 7 in order 
to obtain an approximation for the envelope of T>(t). Over an interval of time of length 691, which is of order e~ 1//2 , 
the value of V at its maximum predicted from equation ( ]27| ) differs from the value obtained by numerical integration 
of the differential equations by only a few percent (actually 2.5%). The error in the comparison of the two values for 
the envelope of V grows slowly over time; in fact, it remains less than 20% over a time interval of length 9000. 

The evolutionary dynamics near resonance is characterized by librational motions described above as well as vari- 
ations in angular momentum. Our numerical experiments suggest considerable variation in the angular momentum 
while the orbit is trapped in resonance. The remainder of this section is devoted to a theoretical explanation of this 
phenomenon. 

During the time that an orbit is captured into resonance, we expect that on average it will reside near an elliptic rest 
point of the first order partially averaged system. Thus, we expect the dynamical variables to satisfy approximately 
equation (^4|) while the orbit is trapped. We will show how to determine the dynamical behavior of the angular 
momentum G under this assumption. In particular, for the (1 : 1) resonance, we claim that G > 0; that is, the orbital 
angular momentum increases. 

It is easy to show that 

dT 1 

-2T s = - — [asm2g - /3cos(2 9 + p)]S m (e), 

dT 1 

-^- + 2T C = — [acos2 5 + /3sin(2 5 + p)]S m (e), 

where 

5 m (e) := m 2 [A m (e) - B m (e)]. 
After substituting these identities as well as the condition (|24|) into the expression for G in (|l9|), we find that 

G = -e[P m {e) + L: 2 S m {e)Q m ], (29) 
where, after some manipulation using the relation e = (1 - G 2 /L 2 ) 1/2 , 

Pm(e) ■= §f [(1 - e 2 ) 3 / 2 (8 + 7e 2 ) 2(8 + ^e 2 + ^e 4 ) 
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Qm '■= - [as'm(mip — 2g) + /3cos(rrup — 2g — p)]. 

The function S m (e) — > as e — ► 0; therefore, it is clear that G > near e = as long as P m (e) < 0. We have 
Fi(e) < -8A/G 7 for < e < 1. For m = 2, we have P 2 (e) < 0, but P 2 (0) = 0, while for m > 2, we have P m (e) > 
near e = 0. Our conclusion is that near the (1 : 1) resonance, we can expect G > during capture provided the 
osculating ellipse is not too eccentric. The other cases are more delicately balanced. However, for a specific choice of 
parameters and for a specific resonance, one can use equation ( p9| ) to predict the behavior of G near the resonance. 

Consider, for instance, the system (Q) with the following parameter values: e = 10~ 4 , A = S/e = 10~ 3 , a = 1, 
j3 = 0, p = 0, and fi = 1. The orbit, as depicted in Figure 0, with initial conditions given by 

(p r ,P9,r,9) = (0.2817,0.6389,1.6628,2.9272) 

appears to be trapped in (1 : 1) resonance with a sojourn time of approximately 40000. Moreover, during this 
period of time the angular momentum appears to increase from approximately 0.66 to 0.95. To demonstrate that 
our theoretical scheme using the second order partially averaged system agrees with our numerical results, let us 
determine the time-scale over which the angular momentum changes from 0.78 to 0.95. According to the numerical 
integration of system (Q) this occurs over about 32700 time units, whereas our approximate analysis of the second 
order partially averaged system predicts a time-scale of 

r 0.95 



1 i-u.ao 

- / [Pi ((1 - G 2 /^) 1 ' 2 )} - 1 dG w 34200, 
e J0.7& 
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a time that is within less than 5% of the value given by our numerical experiment. In this integral, we have used the 
fact that in equation (|29|), the term L~ 2 Si(e)Qi is small compared with Pi(e), 

It should be clear from the results of this section that — by employing the second order partially averaged system — 
we have been able to provide theoretical explanations for the behaviors of semimajor axis a and eccentricity e of 
the orbit that is trapped in resonance. The basic dynamical equations for our model are only valid to order e, since 
physical effects of higher order have been neglected in equation therefore, it would be inappropriate to go beyond 
the second order averaged system in this case. The averaging method is general, however, and can be carried through 
to any order. 



5. CONCLUSION 



Some of the observed structures in the solar system are the direct results of evolutionary dynamics while trapped 
in resonance (see ||, ]l3|, and [jlj]]). We have argued in this paper that such dynamics can in general be explained 
using the method of averaging. In particular, the main aspects of the evolution of the dynamical system during the 
time that it is locked in resonance can be understood on the basis of the second order partially averaged dynamics. 
We have illustrated these ideas using a particular model involving a Keplerian binary system that emits and absorbs 
gravitational radiation. 



APPENDIX A: DELAUNAY EQUATIONS OF MOTION 



The Delaunay equations of motion for our system are given in 0] as 



dL 
~dl 
dG 
~dt 
M 
dt 
dg 
dt 



fdC lf . OS , 
fdC ir . dS ,. , 
/ dC 1 . , dS , , , 



, 'dC , , , 



where 



4>{t) 



-ail 2 cos Qt. 



C and S are given by equation (^|) and 



h 
fa = 
ft = 

f 9 = 



4 
18G 



1 



16 L 



T7 +1 -' 



20. 



9S s 



17 



Sfc, 
+ Sft, 

Sfg, 



1 



pn 2 cos(fif + p), 
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SO 



7G 2 L - - 



1 25 L 2 G 2 



25- 



L 2 G A 



(Al) 



(A2) 



The /d can be expressed, in principle, in terms of Delaunay elements. In fact, expressions of the form r q and 
r~ q sinw can be expanded as Fourier series in terms of the mean anomaly t (see [Q, p. 54). For example, the following 



results are known J? 20 
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L 2 

— = 1 + 2^ J n (ne) cos n£, (A3) 

n=l 

L 4 (l-e 2 ) 1/2 °° 

— ^sinS = y 2nJ n (ne) smn£. (A4) 

n=l 

For the calculation of the other similar expressions — and eventually the fo — see the methods of classical celestial 
mechanics involving the Bessel functions described in JtJ and [^0| . Only the averages of the fr> — calculated in || and 
given in equation (|l^) — are needed for the purposes of this paper. 

APPENDIX B: THE METHOD OF AVERAGING 

In this section, we give a brief introduction to the method of averaging as it is used in this paper. More detailed 
accounts can be found in many excellent sources. We suggest references Q, PH , p2|] , [ jl8| , and jl9[, as they offer a 
range of mathematical sophistication as well as diverse points of view. 

The method of averaging is usually applied to differential equations in one of several standard forms. In celestial 
mechanics, perhaps the most common standard form is "angular standard form" 

i = ef(I,6), 6 = aj(I)+eg(I,6), 

where / is an m-dimensional variable, 9 is an 71-dimensional angular variable, both / and g are 2tt periodic in 9, and 
e is a small parameter. This form is often obtained after a model of disturbed two-body motion is expressed in the 
action-angle variables appropriate for the unperturbed Kepler problem as in equation (M). The second form is "time 
periodic standard form" 

x = eF(x,t), (Bl) 

where x is an m-dimensional variable and F is T-periodic in time. This form is obtained for many different physical 
problems. In particular, in our case — i.e. equation ( p"3| ) — the small parameter is e 1 / 2 and the equation of motion 
is a Taylor series in e 1 / 2 . Of course, the time periodic standard form can be viewed as a special case of the angular 
standard form by introducing a new one dimensional angular variable: 

x = eF*(x,9), 6 = ^. 

However, it is useful to separate the two cases. The Averaging Principle is often applied to differential equations 
after they are transformed to one of the above standard forms. It consists of approximating a solution of the angular 
standard form by the corresponding solution (with the same initial conditions) of the averaged system 

or, for the time periodic standard form, by the equation 

1 f T 

y = e rJ o F ^v,t)dt. 

The validity of the approximation produced by the Averaging Principle is often left unverified. However, as many 
examples show, the Averaging Principle is not valid in general. Thus, a careful analysis of a physical problem must 
always deal with this issue. Fortunately, the problem has been treated in great detail in the mathematical literature. 
Many useful and rigorous results are now available. 

The basic result involves only one angular variable; that is, it pertains either to the time periodic standard form, 
or to the case when 9 is one dimensional in the angular standard form. The Averaging Theorem asserts in this case 
(with appropriate smoothness of the functions involved) that the averaged solution approximates the original solution 
with an error of order e on a time-scale of order 1/e; for example, if x(0) — 1/(0), then there are constants C\ and C2 
that do not depend on e such that 

\x(t) -y(t)\ < Cie, 0<t<^. 
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Furthermore, there is a change of coordinates (the averaging transformation) of the form x = z + eK(z, t) such that 



i r T 

z = e—l F{z,t)dt + e 2 F 2 {z,t) + e s J r (z,t,e). (B2) 
T Jo 

This is exactly the idea that we use to obtain the averaging transformation (|l5|): The purpose of the averaging 
transformation is to render the transformed equation in a form that is already averaged to first order. The second 
order averaged system is then 

ii = e)- [ F(u,t)dt + e 2 ^ [ F 2 (u,t)dt. (B3) 
1 Jo * Jo 



That is, once the differential equation (Bl) is transformed into the exact form (B2) using an averaging transformation 



the second order averaged system is simply obtained from (B2) by averaging the second term and dropping terms 
of higher order. The averaging transformation can be uniquely defined as in Section 2; indeed, this leads to a 
unique second order averaged system. It turns out that a solution u(t,uo,e) of t his system with initial condition 



u(0,uo,e) — uo can be used to approximate the corresponding solution of equation (Bl) with an error of order e on 



a time-scale of order 1/e. More precisely, there are constants c\ > and C2 > such that 

\x(t,u a + eK(u ,0),e) — [u(t,u a ,e) + eK(u(t,u , e),t)}\ < cie 2 

for < t < C2/e. Thus, the second order averaged system gives a second order approximation to the full dynamics. 
We expect that the qualitative features of the dynamics of the original system, for example, the stability of its rest 
points, are reflected in the second order averaged system. While this is not always the case, it is a reasonable working 
hypothesis for the analysis of our model. The partially averaged system is the result of the application of the method 
of averaging to the dynamical system at resonance. Some of the delicate mathematical issues associated with higher 
order averaging, together with an analysis of resonance trapping for a one degree of freedom model, are discussed 
in ||. 

It is important to note that for systems in angular standard form with more than one angle the Averaging Principle 
is not valid in general. In fact, it is exactly the orbits for which the Averaging Principle is invalid that correspond 
to the orbits that are captured into resonances. In ||, we have discussed the implications of the result of Neishtadt, 
which asserts that, under additional assumptions that hold for our model, the set of initial conditions corresponding to 
capture has measure of order e 1 / 2 . Resonance trapping is therefore proved to be rare for the case of two angles unless 
Neishtadt 's hypotheses are violated (see JlTj for a discussion). Rigorous results about the validity of the averaging 
method for the case of more than two angles are much weaker than Neishtadt 's result (cf. p6[). 
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FIG. 1. The plots are for for system (Q) with parameter values e = 10~ 4 , 5/e = 10 -3 , a = 1, /3 = 0, p = 0, and fl = 1. The 
top panel shows L = a 1 ^ 2 versus time for the initial conditions (p r ,pe,r,8) equal to (0.2817, 0.6389, 1.6628, 2.9272). The middle 
panel shows eccentricity versus time and the bottom panel shows G versus L. Here, L* = 1 corresponds to (1 : 1) resonance, 

i/L* = a 



FIG. 2. The plots are for system (Q) with parameter values e = 10 4 , S/e — 10 3 , a = 1, /3 = 0, 
p — 0, and fl = 1. The top panel shows L = a 1 ^ 2 versus time for the initial conditions (p r ,pe,r,6) equal to 
(-0.007485536, 0.850305, 2.758946, 0.97661369282041). The middle panel shows eccentricity versus time and the bottom panel 
shows G versus L. Here, L, w 1.2599 corresponds to (2 : 1) resonance, 2/ Li — fl. 



FIG. 3. The plots are for system (Q) with parameter values e = 10~ 4 , S/e — 10 -3 , a = 1, /3 = 0, p = 0, and fl = 1. The 
top panel shows L = a 1 ^ 2 versus time for the initial conditions (p r ,pe,r, 9) equal to ( — .09951, 1.56904, 2.33813, 2.32909). The 
middle panel shows eccentricity versus time and the bottom panel shows G versus L. Here, L* ~ 1.5874 corresponds to (4 : 1) 
resonance, 4/L 3 = fl. 
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